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Section 1 
INTRODUCTION 


ffc PCM* QUALffY 


1 * 1 Background 

Protein crystals for X-ray diffraction study are usually grown resting on 
the bottom of a hanging drop of a saturated protein solution, with slow 
evaporation to the air in a small enclosed cell. The evaporation rate is 
controlled by hanging the drop above a reservoir of water, with its saturation 
vapor pressure decreased by a low concentration of a passive solute. The drop 
has a lower solute concentration, and its volume shrinks by evaporation until 
the molecular concentrations match. 


Protein crystals can also be grown from a seed crystal suspended or sup- 
ported in the interior of a supersaturated solution. The main analysis of this 
report concerns this case because it is less complicated than hanging drop 
growth. 


Convection effects have been suggested as the reason for the apparent 
cessation of growth at a certain rather small crystal size. It seems that as 
the crystal grows, the number of dislocations increases to a point where fur- 
ther growth is hindered. Growth in the microgravity environment of an orbiting 
space vehicle has been proposed as a method for obtaining larger crystals. 


Experimental 
protein crystals 
others. Figure 1 


observations of convection 
have been reported by Pusey et 
was provided by Pusey. 


effects during the growth of 
al. ( 19Q6 ) , Carter (1987), and 


Figure 1. Observations of 
Convection from a Growing 
Crys ta 1 . 
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1.2 Hydrodynamics of Protein Crystal Growth 


As a protein crystal grows, the protein concentration in the surrounding 
supersaturated solution is decreased. This normally also decreases the den- 
sity, and the fluid rises past the growing crystal and ascends in a buoyant 
p 1 ume . 


This convection flow can be expected to produce the following three 
effects : 


It replaces the surrounding fluid with fresh solution from further 
away. 

Growth is asymmetric, since the concentration depletion is least 
below the crystal, and greatest above it, where the plume leaves. 


The shear flow distribution at the crystal boundary modifies the 
process of crystal growth. 

This third effect is of great concern. In quasi-equilibrium, molecules 
are continually being added to and removed from the crystal, and in the thermal 
agitation, they can settle in a location of minimum energy, forming a quality 
crystal. But if the growth occurs in a shear environment, the relative 
positioning of the added molecules is modified, implying much more dislocations 
in the added layers of crystal. It has been suggested by many workers that the 
accumulation of these dislocations may impede further crystal growth. 


If there is ambient density stratification, due to temperature gradients 
or to concentration gradients produced by the crystal growth or maintained 
externally, the convection flow is modified. In the simplest cases, the am- 
bient supersaturated fluid at a distance from the crystal can be regarded as 
homogeneous . 

There are other hydrodynamic phenomena and issues in hanging-drop protein 
crystal growth. 

One issue is the shape of the hanging drop, and the limits on its 
size. With appropriate support, much larger drops can be used in a 
m icrogravi ty environment. The water reservoir must also be sup- 
ported, of course, for example it can be held by surface tension in a 
porous material. 

A second hydrodynamic issue is the transport of evaporated moisture 
away from the drop, by a combination of diffusion and free or forced 
convect i on. 


The processes of evaporation and condensation lead to density varia- 
tions in the drop and reservoir. The variation of the solute con- 
centration in the reservoir, and of the solute and protein concentra- 
tion in the drop, is controlled by diffusion and convection 
processes , 

A fourth issue is capillary convection in the drop, with surface 
tension variations caused by variations in the temperature and 
concentrat i on. 


These phenomena are discussed further in the appendix. Their detailed 
study is beyond the scope of this present analysis. Attention is confined in 



the rest of this report to the growth of a supported crystal in an large fluid 
vo l ume , 


1.3 Model Assumptions 

We have performed numerical computations of the total flow and concentra- 
tion distributions, including the boundary layer thicknesses and boundary shear 
flow, in the convection near a growing spherical crystal. We made the follow- 
ing assumptions? 

the growing crystal remains spherical; 

its growth rate is negligible compared with the convection flow 
speeds in the surrounding solution; 

the surrounding fluid volume is effectively infinite, and the fluid 
at large distances is homogeneous (except in the ascending buoyant 
p 1 ume ) ; 

temperature variations are negligible; 

effects of the 'sting' support- for the crystal are neglected. 

The limitations of these assumptions are discussed in the following paragraphs. 

The protein crystal will not be spherical. New molecules prefer to add 
themselves to an existing layer, rather than to start a new layer. The crystal 
shapes formed are therefore particular to the individual molecule. In addi- 
tion, the convection flow may establish preferred directions of growth, because 
of variations in flux around the surface and because of velocity shear on the 
sides. This assumption of a spherical shape will of course be better for some 
proteins than for others. 

The growth rate of the crystal, relative to the flow speeds, depends on a 
number of parameters, including the degree of supersaturation of the fluid, the 
protein concentration in the crystal (rarely 100%), and the Rayleigh number 
(discussed below). This assumption is generally a good one. 

The surrounding fluid is effectively infinite if its dimensions are at 
least thirty to a hundred times the crystal radius. The requirements obviously 
depend on the accuracy required. The boundaries modify the flow, and the 
concentration depletion in the plume collects at the top of the chamber and 
spreads downward . 

Temperature effects are negligible because little heat is generated in 
protein crystallization and because the diffusivity for heat is so much more 
than for protein concentration. External temperature effects, and resulting 
convection, may be important in some configurations. 

The support for the crystal will cause minimal flow modification if it is 
a thin ’sting' with the crystal on the end, and if it extends up vertically 
from the bottom of the chamber. A side sting support will cause a larger flow 
modi f icat ion. 
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1.4 Concentration Boundary Condition 


We consider two limiting cases for the concentration boundary condition on 
the crystal surface. In the first limit, the concentration deficiency 
(compared with the supersaturated protein solution at infinity) at the crystal 
surface is fixed, while the growth rate varies with position depending on the 
calculated flux of protein toward the surface. In the second limit, the growth 
rate and flux are fixed, while the concentration deficiency varies over the 
crystal surface. Numerical results are presented for both cases. 

The truth is between these two limits, as illustrated below. The figure 
shows the protein flux (linearly equivalent to the growth rate or to the normal 
derivative of the concentration) as a linear function of the excess concentra- 
tion at the interface over saturation. Note that equilibrium at the interface 
is only a good approximation if growth is extremely slow, because of the large 
size of the protein molecules of interest. 


Growth Rate 



Concentration 
at Interface 


Figure 2. Illustrative Plot of Growth Rate versus Interface Concentration 


In the first limit, the supersaturation at the interface is small compared 
with that at large distances, so that the interface deficiency is effectively 
constant. In the second limit, the supersaturation at the interface is only 
very slightly less than the super saturation at infinity, so that the flux is 
effectively constant. Plainly, the truth will always be somewhere between 
these two extremes. 

We have not included any reduction in the growth rate due to interface 
velocity shear, in our model computations. We have merely computed the shear 
as a function of the radius and the other parameters. It is our hypothesis 
that as the radius increases, the corresponding increase in shear velocity 
modifies the growth process, leading to an increasing number of dislocations. 
Further growth is stopped, since the supersaturation of the fluid cannot be 
increased further without spontaneous nucleation. 
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1.5 Problem Parameters 


The key constant parameters defining the problem are as follows: 
crystal radius a; 
gravity g; 
mean density p; 
kinematic viscosity y5 and 
diffusivity for the protein, D; and either 

surface relative density deficiency C for protein concentration; or 
surface relative density gradient C ? for the second limit. 

Values of interest are as follows: 
a < 0.1 cm, 

2 2 

g = 980 cm/sec (or about 0.01 cm/sec for microgravity), 

P • 1 gm/cc, 

? 

D = 0.01 cm“/sec, 

-7 -5 2 

D = 10 to 10 cnT/sec, and either 
C < 0.0001 (no units), for the first limit, or 
C* < 0.001 cm for the second limit. 


1,6 Dimensionless Parameters 

From these input parameters, we can derive two independent dimensionless 
parameters characterizing the flow. The Schmidt number is 

S = V/D - lo 3 to 10 5 . 

The Rayleigh number is 



for the first limit, and 



for the second. The Grasshof number R/S is also in common use. Using the 
parameter ranges abovg, its maximum possible value is unity, corresponding to R 
values up to about 10 . 
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1.7 Dimensionless Equations 


Using the length scale a and the diffusivity scale D yields the dimension- 
less equations 


S 


-1 


Do 

Dt 





2 

V u_ + Rcz , 


The flow u is zero on the spherical boundary r = 1. The variables 
at infinity. The dimensionless concentration c is unity on r = 1, 
limit where the value is imposed on the growing crystal surface, 
limit of imposing a fixed flux, the boundary condition is 


tend to zero 
in the first 
In the other 


3c/ 3r = -1 . 


To obtain the dimensional variables from the corresponding dimensionless 
variables, they must be multiplied by the following scales: 


1 ength 

speed 

time 

pressure 


a (cm) ; 

D/a ( cm/sec ) ; 

2 

a /D (sec) ; 

2 2 

pD/a“ (dynes/cm ); and 


relative density deficiency C or C'a (no units). 


1.8 Theoretical Estimates 


This system of equations is not amenable to analytic solution, and must be 
solved numerically. The scales of the solution can be estimated, and this is 
an important preliminary to numerical solution. 

In the first place, the momentum term in the equation of motion is always 
small in comparison with the viscous term, since S is so large. 

For R small compared with 1, the concentration solution near the sphere is 
approximated by the conduction (no flow) solution c = 1/r, This applies in 
both limiting cases, with the boundary condition c = 1 or with the boundary 
condition that the radial derivative is -1. The equation of motion determines 
the dimensionless flow speed magnitude as a function of radius, as of order 

R(r-l). 

The corresponding shear on the crystal boundary is of order R, The upper limit 
on the radius, for the conduction solution and the corresponding flow solution 
to apply, is 


The theoretical estimates for the dimensionless flow and boundary layer 
thickness, for large R, are more difficult. Assume boundary layer thickness T 
and velocity shear Q. Then the velocity in the layer is QT. From the concentra- 
tion equation balance, 
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QT = 1/T 2 . 


Also the ratio of the concentration deficiency (compared with infinity) to its 
normal derivative at the interface is T (whichever of these two quantities is 
imposed by the boundary condition). Assuming that the buoyancy in the layer is 
significant in driving the shear (rather than the shear being mostly due to 
plume buoyancy), the motion balance gives in the two limits, 

CIT/T 2 = R, and 

QT/T 2 = RT. 


Thus, when the concentration deficiency as compared with infinity is imposed, 
the shear and thickness are of order 


,3/4 


T = R 


-1/4 


while when the concentration deficiency gradient is imposed, the shear and 
thickness are of order 


Q = 


,3/5 


T 


R 


-1/5 


The numerical results confirm these estimates. 


The plume thickness, as the buoyant plume leaves the crystal, is of the 
same order. The angular thickness decreases without limit as the plume rises 
and its speed increases, so that full resolution of the whole plume using a 
spherical polar mesh is impossible. Fortunately, accurate resolution of the 
far plume and of the far-field flow which it drives is not a requirement, in 
the study of convection effects on protein crystal growth. 

2 

The dimensional shear is QD/a~. Substituting our three formulae for Q, 

Q = R, for R << 1, 

3/4 

Q = R , for R >> 1, fixed deficiency, 

3/5 

Q = R , for R >> 1, fixed growth rate, 


and using our two definitions of R and its dependence on a, gives 

. t u ^ b 

Dimensional shear a , 


where 


b 

b 

b 

b 


1, 

for 

R 

<< 

i, 

2, 

for 

R 

<< 

i, 

1/4, 

for 

R 

>> 

i, 

2/5, 

for 

R 

>> 

i, 


fixed deficiency, 
fixed growth rate, 
fixed deficiency, 
fixed growth rate. 
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More precisely, we can 
shear in the I iraits of sraal 1 
growth rate. 


derive the following table for the dimensional 
and large R, for fixed deficiency and for fixed 


R << 1 
R >> 1 
Table 1. 


Fixed Deficiency 

gCa 

y 

3/4 1/4 

<gC/Ry) (Da) 


Estimates for Dimensional 


Fixed Growth Rate 


gCa“ 

D 

3/5 2/5 

(gC /Dy) (Da) 

Shear Maximum in Four Limiting Cases 
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Section 2 
NUMERICAL METHODS 


2.1 Convection Modelling wit h the AGCE — Code , 


Ue have modeled convection effects in protein crystal growth by using an 
existing computer code. This flexible code was developed by the presen 
author, under MSEC sponsorship, to model thermal convection in a sp • 
layer The design objective for the code was studies related to NASA s 

proposed A+mospheric General Circulation Experiment ( AGCE) . The code was built 
Jith great flexibility, and with minor modifications and appropriately chosen 
input data it was well suited to NASA's requirements in this study. 

The code uses the primitive variables in spherical polar ^ordinates. The 
vorticity-stream function formulation was avoided, to facilitate the use of 

s lity analyses and in three dxmensions. The domain is octangular 
in these coordinates, which allows for the whole region between an inner and 
oute^ sphere? A wide range of options is available for the thermal boundary 
conditions on each boundary segment. With the substitution of concentration 
for temperature, the code was ready for use with almost no other chang. s. 
imposed zero dimensionless density deficiency (scaled concentration deficiency 
on a large outer radius (which is a good approximation as discussed in Section 
1.3). On the inner sphere of radius unity, we imposed one of the two limiting 
dimensionless concentration boundary conditions. 

Th code allows the use of f lex ib 1 y-def ined non-uni form meshes These 

allowed us to resolve the thin boundary layers on the sphere (for Urge , 
litZut wasting mesh points on the regions far from the sphere. Simi ar y, we 
could resolve the thermal plume (at least near the crystal) without jsing 
fine 9 mesh outside the plume. 

The numerical method options available in the AGCE code include iteration 
to a steady state using an algorithm which approximates time-stepping the 
unsteady Rations! but with a time step which varies with position. This was 
a crucial factor in obtaining steady-state solutions in reasonable time, and at 
reasonable cos?in computer resources, since the time scales vary enormously m 
this application over the relevant domain. Implicit algorithms are f ° r 

the advection and diffusion of heat and momentum, and for the representatio 

internal gravity waves. 


2.2 Input Data 
the top line. 

The data reflects the non-dimensional problem, even though the code ac- 
cepts dimensional data. 
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CRYSTAL/S= 1000/R=10000 


•••PROBLEM PARAMETERS* • * 


DOMAIN 
BOUNDARIES 
THL = .OOOOOOO 

THR - 3.1*15926 

RB = 1.0000000 

RT = 30.0000000 

PORR= .0000000 
ZONAL WAVE NUMBER RJ 


BOUNDARY LAYER 
THICKNESS FOR MESH 
DTHL - .0400000 

DTHR = 1.0000000 

DRB = .0800000 

DRT = 30. 1000000 

P0RTH= .0000000 
• 5, 5, 1 


DIFFUSIVITIES 
ROTATION RATES 


ANU 
AKAPPA 
OHEGAC 
OMEGAM 
OMEGA 1 


1000.0000000 
1.0000000 
.OOOOOOO 
.0000000 
.OOOOOOO 


BUOYANCY FORCE 
SURFTC= .OOOOOOO 

GTERR - 1000.0000000 
GDI EL = .OOOOOOO 
PD I EL - .OOOOOOO 
ALPHA = 10000.0000000 
ALPHA2= .OOOOOOO 


INTERIOR HEATING 
POWERS OF TH & R 
A I HEAT = .OOOOOOO 

PIHTTH = .OOOOOOO 

PIHTR = .OOOOOOO 

CIHTTH = .OOOOOOO 

CIHTR = .OOOOOOO 


BOUNDARY TEMPERATURE SIDE & CORNER 

1 - VALUE FIXED TEMPERATURES 

2 - ZERO FLUX * - FIXED, (ALTH.ALR 


3 - AXIS 
LTBOT = 1 
LTTOP = 1 
LTLEFT = 3 
LTRITE = 3 


5 - 
ABOT 
ATOP 
ALEFT 
AR1TE 


FIXED, 
= .00 
= .00 
= .00 
* .00 


> ALTH.ALR 
TTL - .00 
TTR = .00 
TBL * 1.00 
TBR = 1.00 


BOUNDARY TEMPERATURE INTERPOLATION 

LINEAR IN TH» • IT OR R"IT 

IT = -9 St 99 FOR INFINITY 

ALTH = .5000000 ALR » 6.0000000 

ITBOT 1 1 -3:C0S(TH/ALTH) 

ITTOP = 1 -*:L0G(TAN4(TH/2) + T*(A/2>) 
ITLEFT = 1 
I TRITE = 1 


BOUNDARY VELOCITY 

1 - NO - SLIP 

2 - FREE - SLIP 

3 - AXIS 
LUBOT = 1 
LUTOP = 1 
LULEFT = 3 
LURITE = 3 


* - STRESS/SLIP = OM 


OMBOT = .OOOOOOO 
OMTOP = .OOOOOOO 
OMLEFT = .OOOOOOO 
OMRITE * .OOOOOOO 


• ••METHOD PARAMETERS • * • 


NUMBER OF STEPS 
NSTEP = 6127 
INITIAL TEMPERATURE 
TIN1T = .00 


TEMPERATURE TIME STEP 
TAUT - 0.010 

PTTH * 1.000 

PTR * 1.000 


VELOCITY TIME STEP 
TAUU = .00003 
PUTH * 1.000 

PUR * 1.000 


UPUIND DIFFERENCES 
TUPVND -- . 100 

UUPUND = .100 

SMTHU = .000 


REQUIRED CONVERGENCE 
RCONV = 9.000 

DCEVM = 1.000 

SMTHU * .000 


IMPLICIT AMPLITUDES AND CONTROLS 
BETTA - .60 FIXTTH * T 

BETUA = .60 FIXUTH * T 

BETUC = 1.00 FIXUC * T 


IMPLICIT AMPLITUDES AND CONTROLS 
BETTD * .60 FIXTR = T 

BETUD = .60 FIXUR = T 

BETINT * 1.00 FIXINT = T 


PRESSURE METHOD 
NPITER = 7 

PEXTRP =1.00 
ADINT = 0.20 


• ••OUTPUT PARAMETERS'" 


CONCENTRATION C 

DIRECT ACCESS READ A WRITE SEGMENTS 
ISEGR = 1 ISEGW = 1 ISEGS = 1 

(0: ANALYTIC. -1: AFTER PRIOR CASE) — 

I STEP BEGIN A INCREMENT FOR PRINTER DIAGNOSTICS 
VARIABLE PHYSICAL CONTOUR INTEGER CONTOUR PRINT NUMBERS 


WRITE STEP BEGINNING A INCREMENT 
IBEGDA = 0 IINCDA = 900 

(ZERO MEANS NSTEP) RTOP = 

RANGES FOR I AND 


U 

V 

V 
T 
P 

DP 

DU 

DV 

DW 

DT 

PS I 

VORT 

AAH 

OMEG 


9901 9902 
9901 9902 
9901 9902 
6127 9902 
9901 9902 
9901 9902 
8901 9902 
9901 9902 
8901 9902 
8901 9902 
6127 9902 
9901 9902 
9901 9902 
9901 9902 
9901 9902 


9903 9904 
9903 990* 
9903 9904 
6127 9904 
9903 9904 
9903 990* 
6127 9904 
9903 9904 
8127 9904 
6127 9904 
9903 9904 
9903 9904 
9903 9904 
9903 990* 
9903 9904 


9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 
9905 9906 


IB 

1 


IE 

22 

22 

22 

22 

22 

22 

22 

22 

22 

22 

21 

21 

22 

22 

22 


II 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

l 

1 

1 

1 


KB 

22 

22 

22 

22 

22 

22 

22 

22 

no 

LL 

22 

21 

21 

22 

22 

22 


DIAGNOSTIC FREQUENCY CONTOUR LINES 
NDIAG = 0 NCLP * 0 NCLI * 0 
14C0PYS = 1 

PLOTTER IDAOUT BEGIN A INCREMENT 
PHYSICAL INTEGER 
1 9914 9915 9001 

9913 9914 9915 9916 

9913 9914 9915 9916 

1-1 1 9914 9915 9001 

1 9913 9914 9915 9916 

9913 9914 9915 9916 

9913 9914 9915 9916 

9913 9914 9915 9916 

9913 9914 9915 9918 

1 9913 9914 9915 9916 

1 -1 1 9914 9915 9916 

1 -1 9913 9914 9915 9916 

1 -1 9913 9914 9915 9916 

1 -1 9913 9914 9915 9916 

1 -1 9913 9914 9915 9916 


4.00 

K 

KE KI 
1 -1 
1 -1 
1 -1 
1 
1 

1 -1 
l -1 
1 -1 
1 -1 
1 


I NCR 
0 
0 
0 
,0 
0 
,0 
0 
.0 
0 
.0 
0 
.0 
0 
.0 
.0 


Figure 3. 


Input Data for the AGCE Code, 


for a typical case. 
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The left and right e boundaries THL and THR are of course zero and fi, 
corresponding to the upward and downward axes. The inner radius RB is unity. 
The outer radius RB of the computational domain is an approximation to an 
infinite domain? the results should be unaffected by any larger choice. This 
was confirmed using an outer radius of 50. 

The viscosity ANU and the thermal diffusivity AKAPPA are S (1000 in this 
case) and unity. The constants ALPHA and GTERR multiply the temperature to 
give the upward buoyancy acceleration; we therefore always set GTERR to S and 
ALPHA to R. Thus the equation of motion is precisely that in Section 1.7, 
multiplied through by S. 

The code uses nonuniform computational meshes in the radial and theta 
directions. The radial mesh spacing is proportional to the product 

CRT + DRT - r ) <r - RB + DRB ) . 

A similar expression defines the theta mesh, using DTHL and DTHR. The quan- 
tities DRB and DTHL are small, to give adequate resolution of the boundary 
layer on the sphere and of the rising plume. 

The LTxxx variables establish that the "temperature" c is imposed on the 
"bottom" RB and on the "top" RT, while axis boundary conditions are applied on 
the axis ("left" and "right"). The corner temperatures (TTL means "temperature 
top left") determine the imposed temperatures on the top and bottom boundaries 
as zero and unity, independent of the interpolation method. The velocity 
boundary conditions are no-slip (zero) at the top and bottom, and normal axis 
conditions at the axis boundaries. 

The OMxxx quantities are rotation rates, and are of course not applicable 
here. The same applies to other zero problem parameters. 

The first group of method parameters define and control the iteration 
method, with the number NSTEP of iterations, the nominal time steps TAUT and 
TAUU which are used for temperature and velocity at the coarsest part of the 
mesh, and the powers of the mesh spacings in the two directions (PTTH, PTR, 
PUTH, and PUR), used to decrease these nominal time steps. The time steps are 
very small in the plume and in the boundary layer on the inner sphere. The 
upwind differencing parameters introduce a very small amount of upwind dif- 
ferencing in the representations of the temperature and velocity equations, 
only for cases and at mesh points where it is needed, to avoid mesh separation. 
Otherwise accurate central differences are used in the representation of all 
terms. 

The second group of method parameters define and control the implicit 
computation. Implicit methods are required because of the thin layers, high 
speeds, and convergence requirements. The logical FIXxx variables are T 
(true), indicating that implicit methods are used in both coordinate directions 
for both temperature and velocity, and also for internal waves. The BETxx 
variables determine the precise form of the implicit methods used, and values 
between a half and unity are appropriate. 

The first group of output parameters describe the restart options. When 
ISEGR is zero, the computation is begun from analytic initial conditions 
defined by the problem parameters. Complete data for a restart are written to 
a direct access file every I 1 NCDA steps, beginning with segment 1SEGW at step 
IBEGDA (for which zero indicates NSTEP as in this data). With ISEGR equal to 
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1, as in this case, the initial conditions are read in from the previously- 
written segment 1 of the direct access file. 

The remaining output parameters control the printed and plotted output. 
The variable RTGP is the outer radius for plots, and is of course much less 
than the variable RT. Printer graphics and numerical output are produced from 
the run itself, for the indicate variables, beginning at the indicated steps 
(iterations) and at the indicated frequency. For this data, the temperature T 
and the stream function PSI are plotted on the printer at iteration 6127, which 
is the last iteration. The temperature and the changes in the velocity com- 
ponents and temperature are also plotted against the integers, as a measure of 
resolution and convergence. 

The quality plots are produced in a separate computation which uses the 
direct access file for input, and reads this same input data. With this data, 
U, T and PSI are plotted, from the data on segment 1 of the direct access file. 
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Section 3 

FIXED CONCENTRATION DEFICIENCY 


3, 1 Summary of Cases Computed 

Table 2 is a summary of the fixed-deficiency cases which we have computed. 
For S equal to 1000, we did R values increasing by factors of ten from one to a 
million. In addition, we did R values of 1, 1,000, and 100,000, with S equal 
to 100,000. In each case, the table shows two sets of results, obtained using 
a 40x40 mesh and using an 90x80 mesh ( labeled LOW and HIGH RESOLUTION respec- 
tively, to indicate the resolution). For each case and resolution, three 
numerical values are shown: 


the maximum shear Q = 3u ,/ dr , on the interface; and 

t? 

the maximum and minimum normal derivative of the concentration dc/3r, 

In each case, the negative sign is dropped. 

LOW RESOLUTION 
S 

1,000 


HIGH RESOLUTION 


100,000 


R 

SHEAR 

3c/ dr 

max 

dc / dr . 
min 

SHEAR 

dc/ dr 

max 

dc/ dr . 
min 

1 

3.07 

1.491 

0.959 

3.41 

1.542 

0.974 

10 

14.46 

1.997 

0.923 

19.86 

2. 155 

0.922 

100 

66.88 

2.827 

0.860 

70.07 

2.874 

0.870 

1,000 

320.20 

4.413 

0.790 

336.45 

4. 417 

0.810 

10,000 

1539.10 

7.565 

0.703 

1642.20 

7.244 

0.766 

100,000 

7348.10 

13.628 

0.579 

8083.30 

12.726 

0. 714 

1,000,000 

33187.00 

20. 191 

0.216 

39799.00 

23.922 

0.622 

1 

3.07 

1.491 

0.959 

3. 13 

1.501 

0.964 

1,000 

322. 44 

4. 416 

0.791 

332.41 

4.360 

0.811 

100,000 

7609. 10 

13.634 

0,581 

8695.70 

13,868 

0.764 


Table 2. Summary of Numerical Results for the Fixed-Deficiency Cases 


The following points should be noted about these results. 

First, the error of the high-resolution case is less than a quarter of the 
error of the low resolution case, because the method parameters and the domain 
choice RB were improved in each case for the high-resolution run, based on the 
low-resolution results. Comparing the numbers, it can be seen that all the 
high resolution results have satisfactory accuracy. 

Secondly, the differences between the results for the two different values 
of S are very small. For this reason, we did not do an exhaustive analysis of 
the effect of S. This result is consistent with the analysis in Section 1.8, 
since S enters the equations only in the momentum term, which is multiplied by 
the reciprocal of S. 

Thirdly, the minimum of the concentration derivative occurs at the end of 
the plume, and is of relatively little significance. The maximum of the normal 
derivative of the concentration is a measure of the reciprocal of the boundary 
layer thickness, since the interface concentration is unity. 


- 13 - 



Fourthly, the boundary layer thickness becomes very small as R increase, 

- s ?! srtMS - 

increased, in order to resolve the boundary layer and plume. 

Figure 4 is a logarith.Io plot of 

.ra h r;:-'™. r:.:,. .»u » - — «- 

predicted slope in that region of unity. 



Figure 4. Logarithmic Plot 
Function of the Rayleigh Number 


of the Shear Maximum as a 
for the Fixed Deficiency Case. 


tion 

R. 


>ure 5 is a logari th.i o plot of the boundary layer thickness as a fund- 
I « Aeain the slope is close to the predicted value of 1/4 - for 1 * r& ® 
; ; raa n R, the analytic solution T = 1/r is approached, with constant 


thickness unity. 
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LOuCRJ 


Figure 5. Logarithmic Plot of the Thickness Minimum as a Function of 
the Rayleigh Number, for the Fixed Deficiency Case. 


3.2 Plots of Typical Results 

Figures 6 through 9 .re plots of the oonoentr.tior distribution solutions 
o for R values of 1, 100, 10,000 and 1,000,000. The deorease in the ‘Oiokness 
of the boundary layer and plu.e is in.ediately apparent. For s.all R. the 
distortion from the exact no-flow solution c - 1/r is slig . t 

flow is so rapid that the concentration deficiency can only diffuse y 

small distance into the fluid. 

Figures 10 through 13 show the corresponding flow solutions. The stream 

funcUoHrthe volume flux per azimuthal radian, flowing “.“^e fn 
point and the axis. The flow follows the stream lines. The rapid increase in 

the values and in the contour increments should be noted. 

Figures 14 through 17 show the theta component of the velocity, for the 
same cases! The contrast is significant. Resolution problems are apparent at 
the largest Rayleigh number. 
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CONCENTRATION 
UALUE/80/S«l,000/R-100 
MAXIMUM ■ 1.0002 

MINIMUM • 3 . 42961E-04 

INCREMENT « 5.00000E-02 



i 


Figure 7. Concentration for Fixed Def 
with S = 1,000 and R = 100 
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CONCENTRATION 

UALUE/80/S" 1,000/R*1,000,000 
MAXIMUM - 1.0008 

MINIMUM • 7.88059E-04 

INCREMENT - 5.00000E-03 



Figure 9. Concentration for Fixed Deficiency Cas 
with S = 1,000 and R = 1,000,000 


STREAM FUNCTION 
UALUE/80/S- 1 , 000/R- 1 
MAXIMUM - 1.25743E-07 

MINIMUM - -8.2172 

INCREMENT - 0.40000 



Figure 10. Stream Function for Fixed Deficiency Case 
with S = 1,000 and R = 1 
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STREAK FUNCTION 
UALUE/80/S-l,000/R-100 
MAXIMUM • 6 . 18981E-07 

MINIMUM - -72.055 

INCREMENT • 4.0000 



Figure 11. Stream Function for Fixed Deficiency Case 
with S = 1,000 and R = 100 



STREAM FUNCTION 
UALUE/80/S- 1 , 000/R- 10,000 
MAXIMUM - 2.20269E-08 

MINIMUM - -583.44 

INCREMENT - 30.000 


Figure 12. Stream Function for Fixed Deficiency Case 
with S = 1,000 and R = 10,000 




STREAM FUNCTION 
UALUE/80/S- 1,000/R- 1 , 000,000 
MAX I MUM - 2.5662 

MINIMUM « -1253.7 

INCREMENT • 60.000 


Figure 13. Stream Function for Fixed Deficiency Case 
with S = 1,000 and R = 1,000,000 






Fisure 14. Theta Velocity Component for Fi 
with S = 1,000 and R = 1 


Fixed Deficiency Case 


24 - 


SOUTHWARD UELOCITV 
UALUE/80/S -1,000^*100 
MAXIMUM • 2.78777E-02 

minimum • -14.548 

INCREMENT • 0.60000 


Figure 15. Theta Velocity Component for Fixed Deficiency Case 
with S = 1,000 and R = 100 
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SOUTHWARD UEIOCITV 
UALUE/80/S" 1 , 000/R" 1,000,000 


maximum 

minimum 

INCREMENT 


215.47 

-2039.6 

100.00 



Figure 17. Theta Velocity Component for Fixed Deficiency Case 
with S = 1,000 and R = 1,000,000 
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Section 4 
FIXED GROWTH RATE 


4, 1 Summary of Cases 

Table 3 is a summary of the fixed growth rate cases which we have 
computed. For S equal to 1,000, we did R values increasing by factors of ten 
from one to a million. In addition, we did R values of 1, 1,000, and 100,000, 
with S equal to 100,000. In each case, the table shows two sets of results, 
obtained using a 40x40 mesh and using an 80x00 mesh (labeled LOW and HIGH 
RESOLUTION respectively, to indicate the resolution). For each case and 
resolution, three numerical values are shown* 

the maximum shear magnitude; and 

the maximum and minimum concentration. 




LOW 

RESOLUTION 

HIGH 

RESOLUTION 

s 

R 

SHEAR 

c 

max 

c . 
min 

SHEAR 

c 

max 

c . 
min 

1,000 

1 

2.59 

0.890 

0.710 

2. 85 

0.874 

0.690 


10 

10.50 

0.772 

0.542 

13.93 

0.743 

0.508 


100 

40.06 

0.640 

0.409 

41.61 

0.635 

0.403 


1,000 

152.29 

0.500 

0.292 

158.29 

0.496 

0.288 


10,000 

572.76 

0.372 

0. 199 

601.03 

0.368 

0. 197 


100,000 

2121.70 

0.266 

0. 132 

2265. 10 

0.263 

0. 132 


000,000 

7658.60 



8534.50 

0. 182 

0.086 

100,000 

1 

0.26 

0.890 

0.710 


0.874 

0.689 


1,000 

152.89 

0.500 

0.291 

169.42 

0.496 

0.288 


100,000 

2125. 10 

0.274 

0. 136 

2436.80 

0.262 

0. 132 

Table 

3. Summary 

of Numerical 

; Results 

for the 

Fixed Growth 

Rate 1 

Cases 


The following points should be noted about these results. 

First, the error of the high-resolution case is less than a quarter of the 
error of the tow resolution case, because the method parameters and the domain 
choice RB were improved in each case for the high-resolution run, based on the 
low-resolution results. Comparing the numbers, it can be seen that all the 
high reso 1 ut ion results have satisfactory accuracy. 

Secondly, the differences between the results for the two different values 
of S are very small. For this reason, we did not do an exhaustive analysis of 
the effect of S. This result is consistent with the analysis in Section 1: 
since S enters the equations only in the momentum term, which is multiplied by 
the reciprocal of S. 

Thirdly, the maximum of the concentration occurs at the end of the plume, 
and is of relatively little significance. The minimum of the concentration is 
a measure of the reciprocal of the boundary layer thickness, since the normal 
derivative is -1. 

Fourthly, the boundary layer thickness becomes very small as R increase, 
as predicted theoretically in Section 1.8. For this reason, the input data 
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used decreasing values of the parameters DRB and DTHL (cf. Section 2.2) as R 
increased, in order to resolve the boundary layer and plume. 

Figure 18 is a logarithmic plot of the maximum shear on the interface, as 
a function of the Rayleigh number R. The slope at large R is close to the 
predicted value of 3/5. We did not go to small enough R to obtain the 
predicted slope of unity in that region. 


LOQ(SHEAR) 
M r 



LOGCR) 


Figure 18. Logarithmic Plot of the Shear Maximum as 
a Function of the Rayleigh Number, for the Fixed Growth Rate Case. 


Figure 19 is a logarithmic plot of the boundary layer thickness as a 
function of R. Again, the slope is close to the predicted value of 1/5, for 
large R. For small R, the analytic solution T = 1/r is approached, with con- 
stant thickness unity. 
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lou cm 


Figure 19. Logarithmic Plot of the Thickness Minimum as 
a FuncU.no, the B.,leigh Number, (or the Fixed Growth Bate Case. 


A . 2 Plots of Typ ical Results 

IIITvalues lYV'lof, W?Ooi°and°!,MO,Om? e ^he decrease In the 

the boundary layer and plume le immediately apparent. For sm. I B, the d 

‘reoTapir ;sr ™ *„ 

distance into the fluid. 

Figures 24 through 27 show the corresponding flow solutions. The 
function is the volume flux per azimuthal radian, flowing between , he ind .dated 
point and the axis. The flow follows the stream lines. The rapid increase 
the values and in the contour increments should be noted. 


Figures 28 through 31 show the theta component 
same cases. The contrast is significant. Resolution 
the largest Rayleigh number. 


of the velocity, for the 
problems are apparent at 
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CONCENTRATION 
FLUX/80/S- 1 * 000/R- 1 
MAXIMUM • 0.8739-4 

MINIMUM • 2.5-4125E-02 

INCREMENT ■ -4.00000E-02 



Figure 20. Concentration for Fixed Growth Rate Case 
with S = 1,000 and R = 1 
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CONCENTRATION 
FLUX/80/S- 1 , 000/R- 100 
MAXIMUM - 0*63498 

MINIMUM - 2* 13409E-04 

INCREMENT - 3.00000E-02 



Figure 21. Concentration for Fixed Growth Rate Case 
with S = 1,000 and R = 100 



CONCENTRATION 
FLUX/ 80 /S-l*® 0 ®/R - 10 ' 00 ® 
MAXIMUM • 0.36767 

MINIMUM • 3.43907E-05 

INCREMENT - E.00000E-03 


Figure 22. Concentration for Fixed Growth Rate Case 
with S = 1,000 and R = 10,000 




CONCENTRATION 

FLUX/80/S- 1 * 000/R ■ 1 ' 0®® ' ®®® 
MAXIMUM - ® *18177 

MINIMUM • 5.74265E-0S 

INCREMENT ■ 1 »00000E-02 



Figure 23. Concentration for Fixed Growth Rate Case 
with S = 1,000 and R = 1,000,000 
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STREAM FUNCTION 
FLUX/8O/S-1»O00/R - 10® 
MAXIMUM ■ S.1391SE-07 

MINIMUM ■ -49.887 

INCREMENT ■ 3.0000 


Figure 25. Stream Function for Fixed Growth Rate Case 
with S = 1,000 and R = 100 




STREAM FUNCTION 
FLUX/80/S- 1 * 

MAXIMUM - 0.00000E+00 

MINIMUM • -307.54 

INCREMENT ■ 15.000 



Figure 26. Stream Function for Fixed Growth Rate Case 
with S = 1,000 and R = 10,000 
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SOUTHUARD VELOCITY 
FLUX/80/S- 1 , 000/R- 100 
MAXIMUM - 1.88525E-02 

MINIMUM - -9.4702 

INCREMENT - 0.40000 



Figure 29. Theta Velocity Component for Fixed Growth Rate Case 
with S = 1,000 and R = 100 
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SOUTHWARD VELOCITY 
FLUX /80/S" 1>000/R"1 0/000 
MAXIMUM ■ 0. 11398 

MINIMUM - -77.634 

INCREMENT - 3.0000 



Figure 30. Theta Velocity Component for Fixed Growth Rate Case 
with S = 1,000 and R = 10,000 
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SOUTHUARD UELOCITY 
FLUX/80/S- 1 # 000/R- 1 # 000# 000 
NAXIfUin - 7.5812 

niNinun ■ -595.23 

INCREMENT - 30.000 



Figure 31. Theta Velocity Component for Fixed Growth Rate Case 
with S = 1,000 and R = 1,000,000 
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Section 5 
CONCLUSIONS 


Ue have suggested that the main reason for the difficulty in growi g 
protein crystals beyond a certain rather limited size in terrestrial gravity 
the accumulation of crystal defects caused by crystal growth in a shear conve- 
n M would apply both to a hanging drop configuration and to growth 

a supported crystal The configuration of the added molecules is modified 

by the shear flow. 

g. have analyzad the easier case of the growth of a supported protein 

jsr.ur rLainr 8 :;h«!c^r"!tr;nir^tr^t:°::;;uih,ro:^;d 

the flow speeds. We assumed that the surrounding supersaturated protein so u 

Uon is effectively infinite, and ho.ogsosous at infinity, and .. neglected 
sting effects. 

, - -”r\^%e^;nrsh"r:n s re l ;r":;in^ajr:ndi" 

"t" 5 " 

H^h ra ue dei^strated that .ith this non-di.ensional Uation, the 

I'taiSt nunber and the .cento, ter. can be dropped frc. the analysis. 

Ve derived the following table for the di.ensicnal shear in the limits of 
email and large R, for fined deficiency and for fixed grow.h rat . 


R << 1 
R >> 1 


Fixed Deficiency 
gCa 


le C/Ri,> 3 M <ba) 1M 


Fixed Growth Rate 


gCa 

D 

3/5, „ . 2/5 
<gC’ /Di» (Da) 


As a crystal grows, R increases rapidly. !' U 'sMn'ooipa^d'Il' th 

‘“ h r e s hear ^ is°strongly S dependent W on radius, S and therefore Jbe criti- 

U "1 t.d^s U Ielatlve!y insensitive to the other parameters. On the other 
hand, r once R is greater than one, smell changes in the other parameters may 

allow substantial changes in a. 


We solved the equation 
one to a million) and for 
confirmed our analysis. We 
tions over the whole range, 
developed under NASA sponsor 


numerically for Rayleigh numbers in the range from 
jchmidt numbers from 1,000 to 100,000. The resu ts 
produced plots of the concentration and flow solu- 
This work as done using an existing computer code, 
;hip for a different application. 
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appendix 


CONVECTION PROCESSED IN HANGIN G-DROP CRYSTAL GROU IH 


A. 1 Background 


In hanging-drop protein crystal growing 
tal density gradients, and resulting flows, 
in the reservoir fluid. These motions are 
in any one of the three alone would cause 
the other hand, this coupling is weak, and 
tially independent. 


configurations, there are horizon 
in the drop itself, in the air, and 
coupled, so that density gradients 
motion in the other two media. On 
in practice the motions are essen- 


Ax i symmetry is a good approximation, . For r ”dius VanHir 
henispher iqa 1 drop of radius aft), and * •»■>£'*«' ls . z J i: it is 
volu.e !Tb“h. The votuae flu, of we .er rora d , f , dSe5 i n t 0 the "brine* in the 

reservotr , aS cT,« 'TV* -er' " --i/aro ne 8 ,i 8 ibie, so the 
total flux varies only with time. 


A.2 Convection i n the Reservoir 


The water vapor flu. into the brine is Proport iona, ^ to 

the‘ ‘brine aSSU Ther aal ^effect s V *of r IZ condensation are ne 8 'i 8 ib,e, since the 
diffusivity of heat is hi 8 her and the heat can diffuse a y. 

With .ore water vapor arriving near the th^i-eservoir^ This^stab- 

di luted .ore than the brine at the circura er tends to keep the brine con- 

1 ishes a convection flow in the is required to 

centration independent of radius. 0 y water vapor 

do this, particularly as convection in the air tends P 

flux uniform over the brine surface. 


A. 3 Convection in t h e Vapor Phase 

th - air is highest near the drop, and lowest near 
The vapor pressure in the a ^ drQp and heaviest near the brine. 

the brine, so the air g stabilizing Most of the convection is in 

Z S lay^rb^de" iTZoi: 'where the watel vapor spreads horizontally by 

convection, before diffuses dowo to the brine. 


A. 4 Convection in t he Drop, 

The re.oval of the water f f t n J'°esi n and a of S anr secondary solutes 

:rtJ ro pi“ y .ss^ : - 

The protein has a very low < if y ' The concentration variations drive both 

only slowly diffuse through the drop. T hanging drop. The solute- 

surface-tension! ^dr i v^n.^a'har ansoni'c'i rculation in the opposite direction. 
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The following figure shows qua 1 i 
the streamlines of the circulation on t 
drop to the brine causes a low-density 
in the drop and density stratificat 
convection is modified by surface tensi 


ativeiy the isopycnals on the left, and 
e right. The water vapor flux from the 
ay 0 r in the brine, a high-density layer 
on and convection in the air. The drop 
m variations. 



Figure 32, Convection Processes in a Hanging Drop Configuration 
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